# SETWD ===========================
rm(list=ls())


# LIBRARIES ===========================

library(tidyverse)
library(foreign)
library(stargazer)
library(lmtest)
library(multiwayvcov)

# IMPORT DATA ===========================

dfgi <- read.csv("data/sve-08_04_24.csv",stringsAsFactors = F)

# POLARIZATION ===========================

## TAB A11. Logit: Polarization Models ----

temp <- dfgi%>%
  filter(!is.na(v2cacamps))

m1 <- glm(erosion.strict ~ v2cacamps, 
          data=temp, family="binomial")
crse1 <- coeftest(m1, cluster.vcov(m1, temp$country.name))


temp <- dfgi%>%
  filter(!is.na(v2cacamps))

m1 <- glm(erosion.strict ~ v2cacamps + year, 
          data=temp, family="binomial")
crse2 <- coeftest(m1, cluster.vcov(m1, temp$country.name))


temp <- dfgi%>%
  filter(!is.na(v2cacamps)&!is.na(gini_disp))

m1 <- glm(erosion.strict ~ v2cacamps + year + gini_disp, 
          data=temp, family="binomial")
crse3 <- coeftest(m1, cluster.vcov(m1, temp$country.name))


temp <- dfgi%>%
  filter(!is.na(v2smpolsoc)&!is.na(gini_disp))

m1 <- glm(erosion.strict ~ v2smpolsoc + year, 
          data=temp, family="binomial")
crse4 <- coeftest(m1, cluster.vcov(m1, temp$country.name))


temp <- dfgi%>%
  filter(!is.na(v2smpolsoc)&!is.na(gini_disp))

m1 <- glm(erosion.strict ~ v2smpolsoc + year + gini_disp, 
          data=temp, family="binomial")
crse5 <- coeftest(m1, cluster.vcov(m1, temp$country.name))



stargazer(crse1,crse2,crse3,crse4,crse5, title="Logit: Polarization Models",digits=3,
          header=F,df=F,omit.stat=c("f","ser","aic","ll"),
          add.lines = list(c("Observations", nobs(crse1), nobs(crse2), nobs(crse3),
                             nobs(crse4), nobs(crse5))),
          star.cutoffs=c(0.05,0.01,0.001),
          star.char=c("*","**","***"),
          notes="$^{*} p<0.05$; $^{**} p<0.01$; $^{***} p<0.001$",
          covariate.labels=c("Political Polarization",
                             "Societal Polarization","Year",
                             "Gini"),
          notes.append=F,table.placement="h!",
          label="tab:plz",
          out="tables/apptab-plz-10_27_24.tex")

## FIG A3. Inequality and Polarization ----

ggplot(dfgi,aes(x=gini_disp,y=v2cacamps))+
  geom_vline(xintercept = median(dfgi$gini_disp,na.rm=T),
             color="gray50")+
  geom_point(alpha=0.2)+
  stat_smooth(method="loess",color="red")+
  xlab("Gini")+
  ylab("Polarization")+
  theme_bw()

ggsave("figures/gini_polz.png",width=3,height=3)

## FIG A4. Polarization and Democratic Erosion ----

ggplot(dfgi,aes(x=year,y=v2cacamps,group=as.factor(eroder),
                color=as.factor(eroder)))+
  geom_point(alpha=0.2)+
  stat_smooth(method="loess")+
  scale_color_manual(name="Eroder",values=c("#0794cb","red4"),
                     labels=c("No","Yes"))+
  xlab("Year")+
  ylab("Polarization")+
  theme_bw()

ggsave("figures/erod-plz.png",width=5,height=3.5)


## TAB A15. Logit: Polarization and Inequality Interaction ----

temp <- dfgi%>%
  filter(!is.na(gini_disp),!is.na(gdppc),!is.na(prs.ge),!is.na(v2cacamps),
         !is.na(democracy_duration))

m6 <- glm(erosion.strict ~ gini_disp*v2cacamps + log(gdppc) + year + prs.ge + 
            democracy_duration + region, 
          data=temp, family=binomial)
crse1 <- coeftest(m6, cluster.vcov(m6, temp$country.name))
crse1


temp <- dfgi%>%
  filter(!is.na(gini_disp),!is.na(gdppc),!is.na(v2cacamps))

m6 <- glm(erosion.strict ~ gini_disp*v2cacamps + log(gdppc) + year, 
          data=temp, family=binomial)
crse2 <- coeftest(m6, cluster.vcov(m6, temp$country.name))
crse2



temp <- dfgi%>%
  filter(!is.na(gini_disp),!is.na(v2cacamps))

m6 <- glm(erosion.strict ~ gini_disp*v2cacamps, 
          data=temp, family=binomial)
crse3 <- coeftest(m6, cluster.vcov(m6, temp$country.name))
crse3


stargazer(crse3,crse2,crse1, title="Logit: Polarization and Inequality Interaction",digits=3,
          header=F,df=F,omit.stat=c("f","ser","aic","ll"),
          add.lines = list(c("Observations", nobs(crse3),nobs(crse2),nobs(crse1))),
          star.cutoffs=c(0.1,0.05,0.01,0.001),
          star.char=c("\\dagger","*","**","***"),
          notes="$^{*} p<0.05$; $^{**} p<0.01$; $^{***} p<0.001$",
          notes.append=F,table.placement="h!",
          label = "tab:plz-interact",
          out="tables/apptab-plz-interact-11_18_24.tex")


## AUC estimate ----

auc <- function(temp){
  pos.scores <- temp%>%
    filter(!is.na(pred1))%>%
    filter(erosion.strict==1)%>%
    pull(pred1)
  
  neg.scores <- temp%>%
    filter(!is.na(pred1))%>%
    filter(erosion.strict==0)%>%
    pull(pred1)
  
  # estimate AUC
  set.seed(135102)
  reps <- 1000000
  mean(sample(pos.scores,reps,replace=T) > sample(neg.scores,reps,replace=T))
}

temp <- dfgi%>%
  filter(!is.na(gini_disp)&!is.na(gdppc))

m1 <- glm(erosion.strict ~ year + v2cacamps, 
          data=temp, family="binomial")
temp$pred1 <- predict(m1, temp, type="response")
auc(temp)
# AUC: 0.82897

m1 <- glm(erosion.strict ~ year + v2cacamps + gini_disp, 
          data=temp, family="binomial")
temp$pred1 <- predict(m1, temp, type="response")
auc(temp)
# AUC: 0.8773

m1 <- glm(erosion.strict ~ v2cacamps, 
          data=temp, family="binomial")
temp$pred1 <- predict(m1, temp, type="response")
auc(temp)
# AUC: 0.8045

m1 <- glm(erosion.strict ~ v2cacamps + gini_disp, 
          data=temp, family="binomial")
temp$pred1 <- predict(m1, temp, type="response")
auc(temp)
# AUC: 0.8476

m1 <- glm(erosion.strict ~ v2smpolsoc + year, 
          data=temp, family="binomial")
temp$pred1 <- predict(m1, temp, type="response")
auc(temp)
# AUC: 0.8389

m1 <- glm(erosion.strict ~ v2smpolsoc + year + gini_disp, 
          data=temp, family="binomial")
temp$pred1 <- predict(m1, temp, type="response")
auc(temp)
# AUC: 0.8882



